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We have developed a scattering-matrix approach for numerical calculation of resonant states 
and Q-values of a nonideal optical disk cavity of an arbitrary shape and of an arbitrary varying 
refraction index. The developed method has been applied to study the effect of surface roughness 
, and inhomogeneity of the refraction index on Q-values of microdisk cavities for lasing applications. 

' We demonstrate that even small surface roughness (Ar < A/50) can lead to a drastic degradation 

' of high-Q cavity modes by many orders of magnitude. The results of numerical simulation are 

^Nj , analyzed and explained in terms of wave reflection at a curved dielectric interface combined with 

,__( ' the examination of Poincare surfaces of section and Husimi distributions. 
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^ . I. INTRODUCTION 

c/j ' Dielectric and polymeric microcavities represent a great potential for possible applications in lasing optoelecronic 
O ' devicesi*^. In conventional lasers, a significant fraction of optical pump power is lost and a rather high threshold 
> i: power is needed to initiate the lasing effect. In contrast, spherical and disk cavities can be used to support highly 
efficient low-threshold lasing operation. The high efficiency of such devices is related to the existence the natural 
. ,: cavity resonances. These resonances are known as morphology-dependent resonances or whispering gallery modes''. 
The nature of these resonances can be envisioned in a ray optic picture, when light is trapped inside the cavity through 
. ^ the total internal reflection on the cavity-air boundary. 
£^ In dielectric cavities optically pumped quantum wells, wires or dots provide an active medium sustaining the 

^H ^, lasing operation^i^iSiL^. Polymeric microcavity lasers are made with an active medium including host and guest 
O ■ moleculesSiifliii. The absorbed light is transferred from the photoexcited host molecules in the non-radiative way by 
I— I' means of resonant energy transfer to the guest molecules. A stimulated emission from the active medium of dielectric 
and polymeric cavities is trapped in high-Q modes for a very long time. This leads to a significant increase of intensity 

■ of radiation inside the cavity and hence to low-threshold laser operation. 

, One of the most important characteristics of cavity resonances is their quality factor (Q-factor) defined as Q = 

■ 27r* (Stored energy) /(Energy lost per cycle). The high value of the Q— factor results from very low radiative losses 
' that are mainly caused by radiation leakage due to diffraction on the curved interface. An estimation of the Q-factor 

C in an ideal circular disk cavity of a typical diameter d ~ 10/xm for a typical WG resonance gives Q ^ 10^^ (see below, 

' Eq- H^Bfl ). At the same time, experimental measured values reported so far are typically in the range of 10'^ ~ 10'* or 
• lowes^'S'SiS'SiiSiii. A reduction of a Q-factor may be attributed to a variety of reasons including side wall geometrical 
imperfections, inhomogeneity of the diffraction index of the disk, effects of coupling to the substrate or pedestal and 



o 

. others. A detailed study of the effects of the above factors on the characteristics and performance of the microcavity 
. ^ ■ lasers appears to be of crucial importance for the design, tailoring and optimization of Q-values of lasing microdisk 
, cavities. Such the studies would require an effective computational method that can deal with both complex geometry 
^H ^' and variable refraction index in the cavity. 

' One of the most powerful and versatile numerical techniques often used in photonic simulation is the finite difference 
• • . time domain method (FDTD)i2*iii4. A severe disadvantage of this technique in application to the cavities with small 
. 5^ ■ surface imperfections is that the smooth geometry of the cavity has to be mapped into a discrete grid with very small 
lattice constant. This makes the application of this method to the problem at hand rather impractical in terms of 
\^ both computational power and memory. 

5^ ,, Another class of computational methods reduces the Helmholtz equation in the infinite two-dimensional space into 
contour integral equations defined at the cavity boundaries. These methods include the T-matrix technique^^'*^, the 
boundary integral methodsiiiS, and othersiS. These methods are computationally effective and capable to deal with 
the cavities of arbitrary geometry. However, the above methods require the refraction index be constant inside the 
cavity boundary. 

In the present paper we develop a new, computationally effective, and numerically stable approach based on the 
scattering matrix technique that is capable to deal with both arbitrary complex geometry and inhomogeneous refraction 
index inside the cavity. Note that the scattering matrix technique is widely used in analysis of waveguidesSS as well as 
in quantum mechanical simulations^*. This technique was also used for the analysis of resonant cavities for geometries 
when the analytical solution was available^^. 

The main idea of the method consists of dividing the cavity region into N narrow concentric rings. At each i-th 
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boundary between the neighboring rings we calculate the scattering matrix S' that relates the states propagating (or 
decaying) towards the boundary, with those propagating (or decaying) away of the boundary. Successively combining 
the scattering matrixes for all the boundaries^fliSi, S-*- (8) ... (8) S^, we eventually relate the combined matrix to the 
total scattering matrix of the cavity S. In order to calculate the lifetime of the cavity modes (and, therefore their 
Q-factor) we compute the Wigner-Smith lifetime matrij^ which, in turn, is expressed in terms of the total scattering 
matrix SSiSiSi. 

Because at each step we combine only two scattering matrixes, it is not required to keep track of the solution in the 
whole space. This obviously eliminates the need for storing large matrices and facilitates the computational speed. It 
is also well known that the scattering matrix technique (in contrast, for example, to the transfer matrix technique) 
is not plagued by numerical instability, because exponentially growing and decaying evanescent waves are separated 
in course of the computation. Note that the present technique of combining S'-matrixes is conceptually similar to the 
recurrence algorithm for calculating electromagnetic scattering from a multilayered sphere2^*S&. However, in contrast 
to these works, the scattering matrix technique presented here can be applied to the systems where the refraction 
index varies as a function of both radial and angular coordinates. 

The paper is organized as follows. In Section III Al we develop the scattering matrix technique for disk-shaped 
cavities. The results of numerical calculations of resonant states and Q- values of nonideal cavities on the basis of the 
developed technique are presented in Section IIIII We consider and compare two cases, a disk cavity of a constant 
refraction index n with side wall imperfection (surface roughness), and, a disk cavity of an ideal circular shape but 
with inhomogeneous refraction index n — n{r,ip). The results of numerical simulation are analyzed and explained in 
terms of wave reflection at a curved dielectric interface combined with the examination of Poincarc surfaces of section 
and Husimi function. Finally, we present our conclusion in Section ll VI 

II. THE SCATTERING MATRIX APPROACH 
A. Formalism 

We consider a two-dimensional cavity with the refraction index n surrounded by air. Because the majority of 
experiments are performed only with the lowest transverse mode occupied, we neglect the transverse (z-) dependence 
of the field and thus limit ourself to the two-dimensional Helmholtz equation. The two-dimensional Helmholtz equation 
for z-components of electromagnetic field is given by 

{^ + l:§-r + ^^) ^) + ^) = 

where = {Hz) for TM (TE)-modes, and k is the wave vector in vacuum. Remaining components of the 
electromagnetic field can be derived from {Hz) in a standard way. 




FIG. 1: Schematic geometry of a cavity with the refraction index n surrounded by air. The space is divided in three regions. 
In the inner (r < d) and in the outer regions (r > R) the refraction indexes are constant. In the intermediate region d < r < R 
the refraction index n is a function of both r and </?. The intermediate region is divided by A'^ narrow concentric rings. In each 
ring the refraction coefficient is regarded as a function of the angle only, rii — ni{(f). 
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We divide our system in three region, the outer region, r > R, the inner region, r < d, and the intermediate region, 
d < r < R, see Fig. ^ We choose R and d in such a way that in the outer and the inner region the refraction indexes 
are independent of the coordinate, whereas in the intermediate region n is a function of both r and ip. In the outer 
region the solution to the Helmholtz equation can be written in the form 



+00 



(2) 



q— — OQ 



where Hq^\ Hq"^^ are the Hankel functions of the first and second kind of the order q describing respectively incoming 
and outgoing waves. 

We define the scattering matrix S in a standard fashion^ASi , 



B = SA, 



(3) 



where A,B are the column vectors composed of the expansion coefficients Aq,Bq in Eq. The matrix element 

Sqiq — (S)qiq glvcs a probability amplitude of the scattering from the incoming state q into the outgoing state q' . 
Because of the requirement of the flux conservation, the scattering matrix is unitary^, 



SS^ =1, 



(4) 



where I is the identity matrix. The time reversal invariance imposes the symmetry requirement upon the scattering 



matria 



q qq ' 



(5) 



These two conditions can be used to control numerical results for the scattering matrix. 

i-th boundary 



\ A, 



A. / 



i-th strip (i+r)-th strip 

FIG. 2: The intermediate region is divided by A'^ concentric rings of the width 2A; pi is the distance to the middle of the i-th 
ring. States o^a*"*"^ propagate (or decay) towards the i-th boundary, whereas states fe^fc*"*"^ propagate (or decay) away of this 
boundary. The i-th boundary is defined as the boundary between the i-th and (i 4- l)-th rings. 



In order to apply the scattering matrix technique we divide the intermediate region into iV narrow concentric rings, 
see Figs. QEl Within each i-th ring we write down the solution to the Helmholtz equation as a linear superposition of 
the states propagating (or decaying) out of the disk center and the states propagating (or decaying) towards the disk 
center (the detailed form of these states will be given in Section FlIBI see Eq. (|2U|l ). At each i-th boundary (defined 
as a boundary between the i-th and i + 1-th rings) we can introduce the scattering matrix that relates the states 
propagating (or decaying) towards the boundary, {a^} and {a*,^^}, with those propagating (or decaying) away of the 
boundary, {&;„} and {6:^^}, 



,4+1 



1< i < AT- 1, 



(6) 



where a*, 6* are the column vectors composed of the expansion coefficients {a^}, {bl^}, see below, Ea. (|^ . For the 
A''-th boundary between the last iV-th ring and the outer region the scattering matrix S'^ is defined in the form 



B 



= S 



N 



N 



(7) 
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In the inner region {i = 0) the solution to the Helmholtz equation has the form 



+00 



= a°Mnkr)e''^'P, (8) 



q— — oo 



where Jq is the Bessel functions of the order q. For the inner boundary (i — 0) between the inner region and the first 
ring in the intermediate region we define the matrix S° according to 

(9) 



The brief outline of the derivation and the expressions for the scattering matrixes S' are given in Section III CI and 
Appendix ^ 

The essence of the scattering matrix technique is the successive combination of the scattering matrixes in the 
neighboring regions. For example, combining the scattering matrixes for the i-th and i + 1-th boundaries, S' and 
S'"*""*", we obtain the combined scattering matrix 8'''+-'" = S' (g) S^^^ that relates the outgoing and incoming states in 
the rings i and i + aSLSi, 

6-) - ^'^"^(::-) (10) 

— ^11 + ^1281^^ (l — 8228'!^"^ ) S21, 

— S^2 (I ~ 8ll''"S22) S 
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§21^ — ^21 (I 8228]^! ) S21 

^^22 ^ '^22 ' '^21 ''22'^11 ) '^22'^12 

Here and hereafter we use the notation Sn, S12, ... to define the respective matrix elements of the block matrix S. 
Combining step by step all the scattering matrixes for all the boundaries < i < N we numerically obtain the total 
combined matrix 8°^-'^ = S" ® S-*" ® . . . S'^ relating the scattering states in the outer region (i — N) and the states in 
the inner region (z = 0), 

In order to obtain the scattering matrix S defined by Eq. (PJ, we eliminate a from Eq. (|11|1 and find the relation 
between S°'N and S, 



'21 (I 



s^ry's^r+s^r. (12) 



To identify the resonant states of an open cavity we introduce the lifetime matrix (often called as Wigner-Smith 
time-delay matrix)^ 

c ak c ak 

The diagonal elements of this matrix give a time delay experienced by the wave incident in q-th channel and scattered 
into all other channels, 

-l,(fc)-Q..-^E^S,v (14) 
9' 

The delay time t^(A;) experienced by a scattering wave is totally equivalent to the lifetime t — \/2ck of a quasi- 
bound state with complex eigenvector k = k — ik - . It is interesting to note that Smith in his original paper dealing 
with quantum mechanical scatteringS^ chose a letter "Q" to define the lifetime matrix of a quantum system because 
of a close analogy to the definition of a Q-value in electromagnetic theory. The total time delay averaged over all M 
incoming channels can be expressed in the form.'^'^^ 

± f ± 1 T. s) . -L f ^ . , 

9=1 ^ ' IJ,= 1 
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where exp(i6'p) = are the eigenvalues of the scattering matrix S, 9 = X]^=i total phase of the determinant 

of the matrix S, det S = O^lLi -^m ~ exp{i6). 

The resonant states are manifested as peaks in the delay time whose positions determine the resonant wavevectors 
kres, and the heights are related to the Q- value of the cavity according to 

Q = UJToikres)- (16) 



B. Calculation of the wave functions in the intermediate region d < r < R 

In the intermediate region the refraction index n depends on both r and ip. Therefore, in contrast to the inner and 
outer regions, in the intermediate region we can not separate variables and find an exact analytical solution to the 
Helmholtz equation. We can however write down an approximate solution to the Helmholtz equation in each ring. 
For this purpose let us look for the solution in the form ^'(r, ip) = R{r)^{ip). Substituting this solution into Eq. (Q) 
we obtain 

d^R{r) ^ r dR{r) 1 d^^ip) ^ , , 

— k n (r, (p)r . (17) 



R{r) dr"^ R{r) dr ^{f) dip'^ 

Let us now assume that each ring with radius pi has a vanishing width 2 A (see Fig. In this case we can 
regard r as a constant within each i-th ring, r « p^, with the refraction index being a function of the angle only 
n(r, ip) = ni{ip). In this approximation the variables in Eq. H17|l separate such that for i-th ring we can write 



dip'^ 

d^R^Ti) dR'{r. 



{C + kfiip)pl)^\^) = (18) 
CRHr^) = 0, (19) 



drf On 

where C is a constant (which can be both positive and negative), and n = r/pi. The angular function ^''{(p) satisfies 
the cyclic boundary condition $*(0) = $'(27r). The solution of Eq. IjlHI) thus provides an infinite set of eigenvalues 
{Cm} with the corresponding eigenfunctions $^((/3). Generally, Eq. ifTBl) has to be solved numerically. For a given 
eigenvalue the solution of Eq. H19() for the radial wave function can be easily written in the analytical form, and 
the approximate solution to the Helmholtz equation in the i-th ring (situated to the left to i-th boundary) reads 

oo 

^,{n,^) = {<ne(-i+^^^h+l^ei-i-^^^h) $^(^), (20) 

m— 1 



where = {r — pi)/pi and 7^ — \J ^ \ ~ Cm- The states in Eq. H2U|I are grouped according to the convention adopted 

in the previous subsection III Al Namely, the states propagating to the right towards the i-th boundary (e^'^™'"') are 
described by the coefficients whereas the states propagating away from the i-th boundary (e"*'''™''') enter with 

the coefficients {6m}- Note that if 7^ becomes imaginary, 7 = in, the state propagating towards (away of) the z-th 
boundary turns into the states decaying towards (away of) this boundary. 

The wave function ^'^4-1 (r^+i , ip) in the (i + l)-th ring (situated to the right to i-th boundary) is given by the similar 
expression with coefficients Om and bm interchanged, 

00 

*,+i(r,+i,^) ^ ')^~'+^ + a:;rie(-^-'^"^')^'+^) <^^+\^), (21) 

This is because in the {i + l)-th ring the states e+'^"^'^'''+i propagate (or decay) away of the i-th boundary, whereas 
the states e~^"^'^'''+i propagate (or decay) towards the i-th boundary. 



C. The scattering matrix S' at the i-th boundary 

In this section we derive the expression for the scattering matrix S' by matching the wave functions across the 
i-th boundary. Using the condition of the continuity of the tangential components of the electric and magnetic fields 
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at the boundary between two dielectric media, the matching conditions at the i-th boundary (i.e. at the boundary 
between z-th and i + 1 rings ) read 

■^,{r,ip) = *,+i(r,(^), (22) 
1 d^,{r,^) ^ 1 0^,+i(r,(p) 
xfif) dr xl+iif) dr 

where Xiiv) — 1 fo'' TM modes, and Xiif) — k'^^iif) fo^' TE modes. 

In order to derive the expression for the scattering matrix S' in the intermediate region {1 < i < N ^1) we substitute 
the wave functions Eqs. H2UI) . (|21|l into the boundary conditions Eq. (|22|l . Multiplying the obtained equations by 

[^ly^{(p))* and integrating over the angle using the conditions of the orthogonahty J^^ dip ($*„((/?))* ^l^,{ip) — Smm' 
we arrive to two infinite systems of equations for the coefficients aj„, aj^^, bl^, b^^^. After some straightforward algebra 
these systems of equations are reduced to the form prescribed by Eq. ^ with the following result 

S' AKA-^BKA-i. (23) 

The scattering matrixes S*',S'^ (for inner i = and outer i ^ N boundaries respectively) are derived in a similar 
fashion. The expression for S' given by Eq. (|23|l holds for all the boundaries < i < N . A particular form of the 
matrixes A,K, A, B is different for three distinct cases, namely, (a) 0-th boundary (the boundary between the inner 
region (i — 0) and the first ring i = 1 in the intermediate region); (b) i-th boundary, 1 < i < iV — 1, (the boundary 
between i-th and i + 1-th rings in the intermediate region) , and (c) N-th boundary (the boundary between the last 
ring z = A'^ in the intermediate region and the outer region [i = N + 1)). The corresponding expressions for these 
three cases are given in Appendix, Eqs. HA1 |I -HA3 |I . 



III. NONIDEAL MICRODISK CAVITIES 



In this Section we apply the scattering matrix method to the calculation of resonant states and Q- values of nonideal 
microdisk cavities with (a) side wall imperfections and (b) circular cavities with inhomogeneous refraction index 
n = n(r, ip). In order to validate present method, we have also performed numerical calculations for structures where 
the analytical solution was available. This includes, for example, an annular billiard consisting of a dielectric disk 
placed inside a larger disk with some displacement of the disk centev^, as well as an ideal circular disk displaced from 
the origin of coordinate system. In the latter case, the positions of the resonant states and Q-values are obviously 
independent of the choice of the coordinate system. However, from computational point of view this case is not simpler 
than that of a cavity of an arbitrary shape, because the displacement from the origin lifts the radial symmetry and 
makes the separation of variables impossible. As an additional tool to validate the numerical solution we use Eqs. 
Q , jSJ to control the unitarity and symmetry of the scattering matrix. 



A. Ideal circular cavity 

Let us first briefiy analyze the resonant states and Q-values of an ideal circular cavity with the radius R and the 
refraction index n. In this case the scattering matrix can be easily written in analytical form. Employing the matching 
conditions Eq. H22(l between the wave function in the outer region r > R, Eq. and the wave function inside the 
disk given by the Bessel functions Jq {nkr) for r < R, Eq. JH)) , we arrive to the expression for the scattering matrix in 
the formSS 

_ (fer) - e [J:,inkr)/Jq{nkr)] H^^\kr) ^ 

Hq (kr) - e [J'q{nkr)/Jg{nkr)] H\^\kr) 

with ^ = n (l/rt) for TM (TE) modes. Derivatives are taken over the full arguments in the brackets. Resonant states 
of an ideal cavity can be inferred from the scattering matrix Eq. (|24(l using Eq. H15|l . 

Each resonant states of an open disk is characterized by two wave numbers, q and m. These two numbers are directly 
related to the corresponding numbers of the closed resonator of the same radius R. The index m is a radial wave number 
and it is related to the number of nodes of the field components in the radial direction inside the disk. The index 
q is called an angular (or azimuthal) wave number because of the analogy to quantum mechanics where the angular 
momentum is given by Lqm — hq. Equating the quantum and classical angular momenta {Lcias = P^sin^jp — hnk) 
we find the relation between the angular wave number and the angle of incidence x in a classical ray picture^^ 



q = nkR sinx- 



(25) 
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Here we are mostly interested in the whispering gallery modes with high Q- values for which the angle of incidence is 
larger than the angle of total internal reflection, x > Xc (sinxc = l/*^)- For such angles of incidence, the transmission 
probability T of an electromagnetic wave incident on a curved interface of radius p is small, T ^ 1. For the case 
when the radius of curvature is much larger than the wavelength, knp ^ 1, (which applies to majority of cavities), 
the transmission probability reads^i 



T= iTj^lexp 



2 nkp 



3 sin (x) 



(cos^ Xc 



2 \3/2 

cos' X 



(26) 



where Tp is the classical Fresnel transmission coefficient for an electromagnetic wave incident on a flat surface. 




s\nx 

FIG. 3: Transmission coefficient T of a locaffy pfane wave incident on a curved surface with the radii of curvature p as a 
function of the incidence angfe x cafcufated from Eq. 1261 . The angfe of totaf internaf reflection sin^c = 0.56 (corresponding 
to n = 1.8). The inset shows the dependence of the average radius of focaf curvature due to boundary imperfections, p, subject 
to Ar for the present modef of surface roughness. 



Figure 131 iilustrates that T decreases exponentially as the difference x ^ Xc grows. The Q- value of the whispering 
gallery mode g in a cavity of the radius R is related to the transmission probability T, Eq. (|26l) . by the relation^ 



Q 



2nkRcosx 



(27) 



where the classical incidence angle x is related to mode number q by Eq. (|25|l . and T ^ 1. 



B. Nonideal cavities with (a) surface roughness and (b) inhomogeneous refraction index 

In this Section we present the results of numerical calculations of resonant states and Q- values of nonideal cavities. 
We consider separately two cases, (a), a disk cavity of a constant refraction index n but with side wall imperfection 
(surface roughness), and, (b), a disk cavity of an ideal circular shape but with inhomogeneous refraction index 

n = n{r, ip). 

Various studies indicate that a typical size of the side-wall imperfections can vary in the range of 5-300 nm (repre- 
senting a variation of the order of ^--^0.05-1% of the cavity radius) ALSiii. An exact experimental shape of the cavity-air 
interface is however not available. We thus model the interface shape as a superposition of random Gaussian deviations 
from an ideal circle of radius R with a maximal amplitude Ar/2 and a characteristic distance between the deviation 
maxima Al ^ 2t:R/50. In a similar fashion we model the inhomogeneity of the diffraction index in the cavity, where a 
parameter An characterizes a mean deviation of the refraction index n from its average value (n) = 1.8. The variation 
of the refraction index n can be caused by different factors including the presence of quantum well/wires/dots forming 
an active medium of the cavity, the local field intensity dependence n — n{I), and other factors. Examples of typical 
structures under investigation are shown in Fig. ^ 
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FIG. 4: Examples of nonideal cavities studied in the present paper with (a) surface roughness and (b) inhomogeneous refraction 
index, (a) Radius of the disk R = 5^m, n = 1.8, surface roughness Ar = 100 nm. (b) R — 5fim, (n) = 1.8, An — 5%. 

FigureOshows calculated Q-values of the disk resonant cavity for different surface roughnesses Ar and the refraction 
index inhomogeneity An in some representative wavelength interval. Note that we have studied a number of different 
resonances and all of them showed the same trends described below. Besides, here and hereafter we concentrate only 
on TM modes of the cavity, because TE modes exhibit similar features. The calculated dependencies of the Q-values 
on Ar and An are summarized in the insets to Fig. |S| 
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FIG. 5: Dependencies Q = Q(A) for two representative modes TM82,i and TMss^r for the cases of (a) different surface roughness 
Ar and (b) different refraction index inhomogeneities. The values of Ar and An are indicated in Figs, (a) and (b) respectively; 
R = bum, n = 1.8 (a), (n) = 1.8 (b). Note that in the case (b) the resonances shift when An varies. For the sake of clearness 
we plot all the resonances centered around their maxima of the corresponding ideal disk (i. e. An = 0). The broadening of the 
high-Q resonance TMss,? is not discernible on the scale of the figure for all the values of An. Insets in Figs, (a) and (b) show 
the dependencies Q — Q{Ar) and Q = Q(An) respectively. 



Let us first concentrate on the low-Q state TM55 j {q = 55, m = 7). A decrease of the both surface roughness Ar and 
the refractive index inhomogeneity An causes graduate and rather slow decrease of the Q-value of this state as shown 
in the insets to Fig. |S1 This behavior is typical for all other low-Q states. In contrast, the high-Q resonances exhibit 
very different and rather striking behavior. Namely, these resonances show a dramatic decrease of their Q-values even 
for very small values of the surface roughness Ar < A/50. At the same time, the Q-values of the cavity decrease 
much more slowly when the refractive index inhomogeneity An increases. For example, let us choose Ar = 20nm and 
An = 5%. For these values of Ar and An the Q-value of the low-Q state TM55J drops by the same factor of ^ 1.3, 
decreasing from Q « 270 to Q « 205. In contrast, for the very same surface roughness Ar, the Q-value of a high-Q 
state TM82,i drops by the factor of ^ 10^^ decreasing from its value Q « 4 • 10^'^ for an ideal disk to Q ~ 260. At 
the same time, for the above value of An = 5%, the Q-value of this resonance decreases to the value of Q « 1.3 • 10^, 
which corresponds to the drop by the factor ~ 10**. (Note that for the case of an ideal cavity the high-Q resonances 
are so narrow that the numerical resolution does not allow a reliable estimation of their exact values. In this case we 
therefore use Eq. (O to estimate their Q- values.) 
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C. Discussion 



In the previous Section we found that the surface roughness Ar and refraction index inhomogeneity An that produce 
similar degradation of low-Q states, cause strikingly different effect on high-Q resonances. In order to understand these 
features, we shall combine Poincare surface of section and Husimi function methods with an analysis of ray reflection 
at a curved dielectric interface. The Poincare surface of section (SoS) is a powerful tool visualizing the phase space 
for a classical ray dynamics in cavities2iSS. We concentrate on the surface section of the phase space along the cavity 
boundary, r S surf. For a given resonant state with an angular number q, the corresponding ray is launched with the 
angle xo — arcsin (q'/nfci?) according to Eq. if^ . Each reflection at the boundary (characterizing by the polar angle 
if, and the angle of incidence x), corresponds to a single point in the plot. The number of bounces for a given angle 
of incidence xo is chosen in such a way that the total path of the ray does not exceed the one extracted from the 
numerically calculated Q-value for the corresponding resonance, L = ctd = Q/{kn). Figures El (a)-(c), (g)-(i) show 
a Poincare surfaces of section (SoS) for the geometrical rays corresponding to the states TMss^y, TMg2,i for different 
values of the surface roughness Ar. For an ideal circular disk (Ar = 0) the Poincare SoS are obviously straight lines 
corresponding to a constant angle of incidence x — Xo- Figures El (b)-(c), (h)-(i) demonstrate that initially regular 
dynamics of an ideal cavity transforms into a chaotic one even for a cavity with maximum roughness Ar < 20nm. 
In Figs. ini(b)-(c), (h)-(i) Axch approximately indicates the broadening of the phase space due to transition to the 
chaotic dynamics. An important observation is that for a given surface roughness Ar the broadening of the phase 
space is independent of angular mode q, i.e. it is the same for low- and high-Q states. 
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FIG. 6: Poincare surfaces of section for geometrical rays corresponding to the states g = 55 (a)-(c) and q = 82 (g)-(i) for the 
cavity with the surface roughness Ar = 0, 20nm, lOOnm. The Husimi distributions for the states TMss^y (d)-(f) and TMg2,i 
(j)-(l) for the same values of Ar as in the corresponding Poincare SoS. Axch indicate the broadening of the phase space due to 
transition to the chaotic dynamics. Dashed lines show the angle of total internal reflection Xc 



We complement classical Poincare SoS by the Husimi function analysis2i22i22.. The Husimi function (often called 
also Husimi distributions) H{ip, x) represents a quantum (wave) analog to a classical Poincare SoS. It is defined as a 
projection of a given cavity mode ^'(r e surf, tp) taken at the surface of cavity into a Gaussian wave packet <&(</?'; ^, x) 
impinging the cavity boundary with the coordinate tp at the angle x, 

H{^, X) = / d^'^ir e surf, ^'M^'; ^, x), (28) 
Jo 
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where the minimum-uncertainty wave packet centered around ip, x with the dispersion in position \J a jl is given by 

1 



2a 



{(p' ~ ip + 2ttI) - ik sin x{(p + ^irl) 



(29) 



where we have chosen a = \/2/k. The Husimi distributions, Figs. ini(d)-(f), (j)-(l), exhibit the same trends as the 
classical Poincare SoS. Indeed, broadening the phase space with the increasing the surface roughness Ar for the 
Husimi functions is the same as the corresponding broadening A^ch in the Poincare SoS. (Illustrative examples of the 
wave functions in cavities for different surface roughness An are shown in Fig. d). 




FIG. 7: Illustrative examples of intensity distribution for the resonant state TMss.y in cavities with Ar = (a), Ar = 20nm 
(b), Ar — 20nm (c); R = n — 1.8. Dashed lines indicate boundaries of the cavity. 




FIG. 8: The Husimi distributions for the states TMss^y (a) and TMg2,i (b) for the cavity with the refraction index inhomogeneity 
An = 5%. 

Figure |S1 shows the Husimi distributions for a circular cavity with an inhomogeneous refraction index. The variation 
of the refraction index An = 5% is chosen in such a way that the degradation of the Q- value for the low-Q resonance 
TM55J is the same as the one for the case of surface roughness Ar — 20nm shown in Fig. As expected, the 
broadening of the Husimi distribution due to increase of An is of the same order as for the corresponding values Ar 
(compare Figs. |Sland|Sl). 

According to Eq @, one can expect an increase of the transmission coefficient (and therefore decrease of the Q- value 
of the cavity) due to the broadening of the phase space Axch, because the incidence angle x effectively moves closer 
to the angle of the total internal reflection Xc- ATch in Fig. 01 indicates the estimated increase of the transmission 
coefficient due to the broadening of the phase space, Axch, as extracted from the Poincare SoS for the Ar = 20nm 
and An = 5%. For the low-Q resonance TMss^y this corresponds to the decrease of the Q-value by the factor of 
AQch ^2^ch^ ~ 1-^' which is consistent with the calculated decrease of the low-Q resonances. 

For the case of high-Q resonance TM82.1 the estimated decrease of the Q-factor is AQch ^^ch^ ~ (^^^ 
Fig. 121), which is consistent with the calculated decrease of this resonance for the case of the inhomogeneous refraction 
index only. (Note that because of a rather approximate definition of A^ch we can give only very rough estimation 
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of the factor ATch-) On the contrary, for the case of high-Q resonances in the presence of surface imperfections, this 
estimated value of AQch is in many orders of magnitude smaUer that the actual calculated decrease of the Q-factor 
(given by factor of w 10^^, see Fig. EJ. 

To explain the rapid degradation of high-Q resonances, we concentrate on another aspect of the wave dynamics. 
Namely, the imperfections at the surface boundary effectively introduce a local radius of a surface curvature p that 
is distinct from the disk radius R (see illustration in Fig. 0)). One may thus expect that with the presence of the 
local surface curvature, the total transmission coefficient will be determined by the averaged value of p rather than 
by the disk radius R. The dependence of p on surface roughness Ar for the present model of surface imperfections 
is shown in the inset to Fig. |3| Figure |31 demonstrates that the reduction of the local radius of curvature from 5/xm 
(ideal disk) to 1.7/zm (Ar = 20nm) causes an increase of the transmission coefficient by ATcur ~ 10^. This estimate, 
combined with the estimate based on the change of ATch is fully consistent with the actual computed decrease of the 
Q-factor shown in Fig. [S] We thus conclude that the main mechanism responsible for the rapid degradation of high-Q 
resonances in non-ideal cavities is the enhanced radiative decay through the curved surface because the effective local 
radius (given by the surface roughness) is smaller that the disk radius R. 

In contrast, the degradation of low-Q resonances (as well as high-Q resonances in the case of inhomogeneous 
refraction index only), is mostly related to the broadening of the phase space caused by the transition to the chaotic 
dynamics. It should be noted however that both factors (broadening of the phase space and the enhancement of the 
transmission due to decrease of the effective radius of curvature) may play a comparable role in degradation of the 
low-Q whispering-gallery resonances in the presence of surface roughness. 

It is interesting to note that an analogues degradation of high-Q modes was recently found in hexagonal-shaped 
microcavities, where the modes were strongly influenced by roundings of the corners even when the characteristic length 
scale (the local radius of curvature) was one order of magnitude smaller than the wavelengthSi. It is worth mentioning 
that one often assumes that long-lived high-Q resonances in idealized cavities (e.g. in ideal disks, hexagons, etc.) are 
not important for potential application in optical communication or laser devicesiSiS^ because of their extremely 
narrow width. Our simulations demonstrate that it is not the case, because in real structures the Q-values of these 
resonances becomes comparable to those of intermediate-Q resonances already for small or moderate surface roughness 
of Ar ~ 10 - 50 nm. 



IV. CONCLUSIONS 

In the present paper we develop a new, computationally effective, and numerically stable approach based on the 
scattering matrix (^-matrix) technique that is capable to deal with both arbitrary complex geometry and inhomo- 
geneous refraction index inside the two-dimensional cavity. The derivation is based on the separation of the cavity 
region into N narrow concentric rings and calculation of the ^-matrix at every boundary between the rings. The total 
5- matrix is obtained in a recursive way by successive combination of the scattering matrixes for all the boundaries. 
In order to calculate the lifetime of the cavity modes (and, therefore their Q-factors) we compute the Wigner-Smith 
time delay-matrix which, in turn, is expressed in terms of the total scattering matrix. 

We apply the developed algorithm to the calculation of resonant states and Q-values of nonideal microdisk cavities 
with (a) side wall imperfections and (b) circular cavities with inhomogeneous refraction index n = n{r,Lp). We find 
that the surface roughness Ar and refraction index inhomogeneity An that produce similar degradation of low-Q 
states, cause strikingly different effect on high-Q resonances. In particularly, in the case of inhomogeneous refraction 
index the increase of An causes rather graduate decrease of the Q-value of high-Q resonances. In contrast, in the 
presence of surface roughness even small imperfections (Ar < A/50) can lead to a drastic degradation of high-Q cavity 
modes by many orders of magnitude. 

In order to understand these features, we combine Poincare surface of section and Husimi function methods with an 
analysis of ray reflection at a curved dielectric interface. We argue that the main mechanism responsible for the rapid 
degradation of high-Q resonances in non-ideal cavities with the surface roughness is the enhanced radiative decay 
through the curved surface because the effective local radius (given by the surface roughness) is smaller that the disk 
radius R. In contrast, the degradation of low-Q resonances (as well as high-Q resonances in the case of inhomogeneous 
refraction index only), is mostly related to the broadening of the phase space caused by the transition to the chaotic 
dynamics. 
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APPENDIX A: EXPRESSIONS FOR THE MATRIXES S' 



In this appendix we present the expressions for the matrixes A, K, A, B entering Eq. (|23|l for the scattering matrix 
S' relating incoming and outgoing states at the i-th boundary. We distinguish three different cases as specified below. 
(a) 0-th boundary (the boundary between the inner region (i = 0) and the first ring i = 1 in the intermediate region), 



All=I, iA22),nj=e-i^'d,nj, A12 = A21 = 0, (Al) 

Kii = I, (K22)™j = e''''"^'Sr,,j, K12 = K21 = 0, 
A = TTfi 1^1 , B 



(J)mj = Jminokd) 6,nj, (J')mj = Jnii'^okd) Smj 



(V 



nakpi J a xim 



(b) i-th boundary, 1 < i < N — I. (the boundary between i-th and i + l-th rings in the intermediate region) 

(Aii)„„ - e^^*(5™j, (A22)mj =e-^^'+i(5„j, A12 = A21 = 0, (K)^, ^ e'^^^^Smj (A2) 



Qi ui,i+ipi+i y ' \ypi -u'^'+iQ'+i 

('cj A'^-t/i boundary (the boundary between the last ring i — N in the intermediate region and the outer region 
(i = N+l)) 



(Aii)„.j ^ e^^"(5„y, A22 = I, A12 = A21 = (A3) 
(Kii)„, = e'^™^",5™j, K22 = I, Ki2 = K2i = 0, 



X%+i{^) 



In Eqs. ljAl|l - ljA3|l the matrixes Q',P' are defined according to 



(P")™, = ( - o + ) -5™,, (Q')™, - ( - ^7™ ) 1 < * < A^. 

J^, Hm^\ and J^', Hm'^^ are the Bessel and Hankcl functions and their derivatives, and — A/pi. 
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